home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / scheme / scm / jackal1a.lha / jacal / norm.scm < prev    next >
Encoding:
Text File  |  1992-12-24  |  7.0 KB  |  224 lines
  1. ;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
  2. ;;; Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.
  3. ;;; See the file "COPYING" for terms applying to this program.
  4.  
  5. (proclaim '(optimize (speed 3) (compilation-speed 0)))
  6.  
  7. (define (vsubst new old e)
  8.   (cond ((eq? new old) e)
  9.     ((number? e) e)
  10.     ((bunch? e) (map (lambda (e) (vsubst new old e)) e))
  11.     ((var_> new (car e)) (univ_norm0 new (cdr (promote old e))))
  12.     ((var_> old new)
  13.      (poly_resultant (list old (list new 0 -1) 1) e old))
  14.     (else (poly_resultant (list new (list old 0 -1) 1) e old))))
  15.  
  16. ;;; Makes an expression whose value is the variable VAR in the equation
  17. ;;; E or (if E is an expression) E=0
  18. (define (suchthat var e)
  19.   (set! e (poly_subst0 _@ (licit->poleqn e)))
  20.   (extize (normalize (vsubst _@ var e))))
  21.  
  22. ;; canonicalizers
  23. (define (normalize x)
  24.   (cond ((bunch? x) (map normalize x))
  25.     ((symbol? x) (eval-error "normalize symbol? " x))
  26.     ((eqn? x)
  27.      (poly->eqn (signcan (poly_square-and-num-cont-free
  28.                   (alg_simplify (eqn->poly x))))))
  29.     (else (expr_normalize x))))
  30. (define (expr_normalize p)
  31.   (if (expl? p) (set! p (expl->impl p)))
  32.   (expr_norm-or-signcan
  33.    (alg_simplify (poly_square-free-var
  34.           (if (rat? p) (alg_clear-denoms p) p)
  35.           _@))))
  36. (define (extize p)
  37.   (cond ((bunch? p) (eval-error "cannot suchthat a vector" p))
  38.     ((eqn? p) p)
  39.     ((expl? p) p)
  40.     ((rat? p) p)
  41.     (else
  42.      (set! newextstr (sect:next-string newextstr))
  43.      (let ((v (defext (string->var
  44.                (if (clambda? p) (string-append "@" newextstr)
  45.                    newextstr))
  46.             p)))
  47.        (set! var-news (cons v var-news))
  48.        (var->expl v)))))
  49.  
  50. ;; differentials
  51.  
  52. (define (total-diffn p vars)
  53.   (if (null? vars) 0
  54.     (poly_+ (poly_* (var->expl (var_differential (car vars)))
  55.             (poly_diff p (car vars)))
  56.         (total-diffn p (cdr vars)))))
  57. (define (total-chain-exts drule es)
  58.   (if (null? es) drule
  59.       (total-chain-exts (poly_resultant
  60.              drule
  61.              (total-diffn (extrule (car es))
  62.                      (poly_vars (extrule (car es))))
  63.              (var_differential (car es)))
  64.             (union (cdr es) (alg_exts (extrule (car es)))))))
  65. (define (total-differential a)
  66.   (cond
  67.    ((bunch? a) (map total-differential a))
  68.    ((eqn? a) (poly->eqn
  69.           (total-diffn (eqn->poly a) (poly_vars (eqn->poly a)))))
  70.    (else
  71.     (let ((aes (alg_exts a)))
  72.       (if (and (null? aes) (expl? a))
  73.       (total-diffn a (poly_vars a))
  74.       (let ((pa (licit->poleqn a)))
  75.         (total-chain-exts
  76.          (vsubst _@ d@ (poly_resultant
  77.                 pa (total-diffn pa (poly_vars pa)) _@))
  78.          aes)))))))
  79.  
  80. (define (diff a var)
  81.   (cond
  82.    ((bunch? a) (map (lambda (x) (diff x var)) a))
  83.    ((eqn? a) (poly->eqn (diff (eqn->poly a) var)))
  84.    (else
  85.     (let ((aes (alg_exts a)))
  86.       (if (and (null? aes) (expl? a))
  87.       (poly_diff a var)
  88.       (let* ((td (total-differential a))
  89.          (td1 (app* _@1/@2 td (var->expl (var_differential var))))
  90.          (vd (var_differential var))
  91.          (dpvs '()))
  92.         (poly_for-each-var
  93.          (lambda (v)
  94.            (if (and (not (eq? vd v))
  95.             (var_differential? v))
  96.            (set! dpvs (adjoin v dpvs))))
  97.          td)
  98.         (reduce-init (lambda (e x) (poly_coeff e x 0)) td1 dpvs)))))))
  99.  
  100. ;;;; FINITE DIFFERENCES
  101. ;;; shift needs to go through extensions; which will create new
  102. ;;; extensions (yucc).    It is clear what to do for radicals, but other
  103. ;;; extensions will be hard to link up.  For instance y: {x|x^5+a+b+9+x}
  104. ;;; needs to yield the same number whether a or b is substituted first.
  105. (define (shift p var)
  106.   (vsubst var
  107.       _@2
  108.       (poly_resultant (list _@2 (list var -1 -1) 1)
  109.               p
  110.               var)))
  111. (define (unsum p var)
  112.   (app* _@1-@2 p (shift p (explicit->var var))))
  113. (define (poly_fdiffn p vars)
  114.   (if (null? vars) 0
  115.     (poly_+ (poly_* (var->expl (var_finite-differential (car vars)))
  116.             (unsum p (car vars)))
  117.         (poly_fdiffn p (cdr vars)))))
  118. (define (total-finite-differential e)
  119.   (if (bunch? e)
  120.       (map total-finite-differential e)
  121.     (poly_fdiffn e (alg_vars e))))
  122.  
  123. ;;;logical operations on licits
  124. ;(define (impl_not p)
  125. ;  (poly_+ (poly_* (licit->poleqn p)
  126. ;          (var->expl (sexp->var (new-symbol "~")))) -1))
  127.  
  128. ;(define (impl_and p . qs)
  129. ;  (cond ((bunch? p) (impl_and (append p qs)))))
  130.  
  131. (define (expl_t? e) (equal? e expl_t))
  132. (define (ncexpt a pow)
  133.   (cond ((not (or (integer? pow) (expl_t? pow)))
  134.      (deferop '^^ a pow))
  135.     ((eqns? a) (math-error "Expt of equation?: " a))
  136.     ((not (bunch? a)) (fcexpt a pow))
  137.     ((expl_t? pow) (transpose a))
  138.     (else (mtrx_expt a pow))))
  139.  
  140. ;;;; Routines for factoring
  141. (define (poly_diff-coefs el n)
  142.   (if (null? el)
  143.       el
  144.     (cons (poly_* n (car el))
  145.       (poly_diff-coefs (cdr el) (+ 1 n)))))
  146. (define (poly_diff p var)
  147.   (cond ((number? p) 0)
  148.     ((eq? (car p) var) (univ_norm0 var (poly_diff-coefs (cddr p) 1)))
  149.     ((var_> var (car p)) 0)
  150.     (else (univ_norm0 (car p) (map-no-end-0s
  151.                    (lambda (x) (poly_diff x var))
  152.                    (cdr p))))))
  153. (define (poly_diff-all p)
  154.   (let ((ans 0))
  155.     (do ((vars (poly_vars p) (cdr vars)))
  156.     ((null? vars) ans)
  157.     (set! ans (poly_+ (poly_diff p (car vars)) ans)))))
  158.  
  159. ;;;functions involved with square free.
  160. (define (poly_square-free-var p var)
  161.   (poly_/ p (poly_gcd p (poly_diff p var))))
  162.  
  163. (define (poly_square-and-num-cont-free p)
  164.   (if (number? p) (if (zero? p) p 1)
  165.       (poly_* (poly_square-and-num-cont-free (univ_cont p))
  166.           (poly_/ p (poly_gcd p (poly_diff p (car p)))))))
  167.  
  168. (define (poly_factorq p) (poly_factor-all p))
  169.  
  170. (define (poly_factor-var c var)
  171.   (poly_factor-split c (poly_diff c var)))
  172.  
  173. (define (poly_factor-all c)
  174.   (poly_factor-split c (poly_diff-all c)))
  175.  
  176. ;;; This algorithm is due to:
  177. ;;; Multivariate Polynomial Factorization
  178. ;;; David R. Musser
  179. ;;; Journal of the Association for Computing Machinery
  180. ;;; Vol 22, No. 2, April 1975
  181.  
  182. (define nreverse reverse)        ;needs to be defined
  183.  
  184. (define (poly_factor-split c splitter)
  185.   (let ((d '()) (aj '()) (b (poly_gcd c splitter))) ; changed #f's to ()
  186.     (do ((b b (poly_/ b d))
  187.          (a (poly_/ c b) d))
  188.         ((number? b)
  189.          (if (eqv? 1 b)
  190.              (cons a aj)                ; nreverse removed
  191.              (cons b (cons a aj))))     ; nreverse removed
  192.       (set! d (poly_gcd a b))
  193.       (set! a (poly_/ a d))             ; 'a' is used as a temporary variable
  194.       (if (not (eqv? a 1))              ; keeps extraneous `1's
  195.         (set! aj (cons a aj))))))       ;  out of the results list
  196.  
  197. ;;; the following algorithm attempts to separate factors in a multivariate
  198. ;;; polynomial with major variable.  It substitues 0 for each variable
  199. ;;; that it finds in turn and takes GCD against the original expression.
  200. ;;; It assumes that it's argument is squarefree and contentfree in the
  201. ;;; major variable.
  202. (define (univ_split pe varlist)
  203.   (cond ((unit? pe) (list))
  204.     ((null? varlist) (list pe))
  205.     ((let ((p0 (signcan
  206.             (poly_gcd pe (poly_subst0 (car varlist) pe))))
  207.            (cvl (cdr varlist)))
  208.        (if (unit? p0)
  209.            (univ_split pe cvl)
  210.          (nconc (univ_split (poly_/ pe p0) cvl)
  211.             (univ_split p0 cvl)))))))
  212.  
  213. (define (univ_split-all poly) (univ_split poly (poly_vars poly)))
  214.  
  215. (define (factor-free-var p var)
  216.   (poly_gcd p (poly_subst0 var p)))
  217.  
  218. (define (factor_test)
  219.   (test (list (list x -1 -2) (list x -1 0 1))
  220.         poly_factorq
  221.         (list x -1 -4 -3 4 4)))
  222.  
  223. ;;;    Copyright 1989, 1990, 1991, 1992 Aubrey Jaffer.
  224.